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ABSTRACT 

We consider the dynamical evolution of two planets orbiting in the vicinity of a first order mean motion reso¬ 
nance while simultaneously undergoing eccentricity damping and convergent migration. Following Goldreich 
& Schlichting|lf2()14), we include a coupling between the dissipative semimajor axis evolution and the damping 
of the eccentricities. In agreement with past studies, we find that this coupling can lead to overstability of the 
resonance and that for a certain range of parameters capture into resonance is only temporary. Using a more 
general model, we show that whether overstable motion can occur depends in a characteristic way on the mass 
ratio between the two planets as well as their relative eccentricity damping timescales. Moreover, we show that 
even when escape from resonance does occur, the timescale for escape is long enough such at any given time a 
pair of planets is more likely to be found in a resonance rather than migrating between them. Thus, we argue 
that overstability of resonances cannot singlehandedly reconcile convergent migration with the observed lack 
of Kepler planet pairs found near resonances. However, it is possible that overstable motion in combination 
with other effects such as large scale orbital instability could produce the observed period ratio distribution. 

Subject headings: celestial mechanics - planets and satellites: dynamical evolution and stability 


1. INTRODUCTION 

Analyses of the Kepler data which take into account obser¬ 
vational and instrumental biases indicate that Sun-like stars 
and stars less massive than the Sun commonly host planets 
(Fressin et akj|2013[ jDressing & Charbonneau||2013| |Mor- 
ton & Swift||2014|. The known planets tend to be smaller 
than Neptune, orbit with periods less than ~ 100 days, and 
a significant fraction of them reside in multi-planet systems 


(Batal hTeUd.|2013t|Munally et al.|2015)|Rowe et al.|2015 1 . 
Given that these planetary systems may represent the domi¬ 
nant outcome of planet formation, advancing our understand¬ 
ing of their formation and past dynamical evolution is a major 
goal of exoplanet science. 

Although these planets typically have weakly constrained 
masses, densities, and orbital elements, there are several 
broad features of the population which can provide clues as to 
the formation of these systems. One interesting characteristic 
is that the period ratios of pairs of adjacent planets do not pref¬ 
erentially lie near mean motion resonances. Initially, this was 
taken as a clue that large-scale migration caused by dissipa¬ 
tive interaction with a gaseous protoplanetary disk (e.g. Kley 
|& Nelson|2012| |Baruteau et al.|20i~4| ) did not act in a particu¬ 
larly important way, since rudimentary models of convergent 
migration between planets robustly predict capture into res¬ 
onance assuming the migration rate is slow enough and the 
eccentricities are small enough. This apparent contradiction 
between migration models and the observed data led to the 
idea of in-situ formation, where these systems ultimately as¬ 
sembled via the same mechanisms that created the terrestrial 
planets in our Solar Syste m, but ac tin g at sig nificantly smaller 
orbital distances (e.g. Hansen & Murray|2012| |2013( Chiang 
& Laughlin|20l3). 

A second feature of these planets complicates this narra¬ 
tive. By modeling the composition of the subset of planets 
with measured masses and radii, Rogers ( }2015j > showed that 
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the majority of planets with radii larger than ~ 1.6f?Q require 
significant gaseous atmospheres. Processes which produce 
volatiles after the formation of a planet (such as out-gassing) 
are thought to proceed at a slower rate than processes which 
strip a planet of its volatile s (like photo-evaporation) ( |Lopez| 
|et al.l|20T2| |Rogers||2015| l. This suggests that the fraction 
of volatile rich planets was larger in the past compared with 
the observed sample, and it also indicates that these planets 
formed while the gaseous protoplanetary disk was present. 

Therefore, in-situ formation no longer provides an immedi¬ 
ate explanation for the observed lack of pairs near mean mo¬ 
tion resonances by simply obviating the need for large-scale 
migration. That is, irrespective of where close-in small plan¬ 
ets originated, they would have interacted with their gas-rich 
natal disks. Therefore it is necessary to find a way to explain 
the lack of observed near commensurabilities in the context 
of planetary evolution within a gaseous disk, the presence of 
which seems required to explain the volatile rich nature of a 
subset of the planets observed. 

Both turbulence in the disk and s mall residual eccentr icities 
can prevent capture into resonance (Adams et al.|2008 Rein 


2012| [Paardekooper et al.|2013| |Batygin|2015)>. Additionally, 

recent work has suggested that a more complete treatment of 
convergent migration and eccentricity damping alone can ac¬ 


count for the lack of pairs near resonance. Particularly, Gol¬ 


dreich & Schlichting (2014) show that in this more complete 


model, which takes into account how eccentricity damping af¬ 
fects the semimajor axis evolution, the long-term stability of 


resonances can be compromised (see also Meyer & Wisdom 
2008). That is, capture into resonance occurs but it is only 
temporary, and escape from resonance occurs on timescales 
comparable to the eccentricity damping time. Since this is 
small compared with the time spent migrating between res¬ 
onances (the semimajor axis decay timescale), the expected 
result is a distribution of period ratios which disfavors reso¬ 
nant valued 


3 These results were derived for the circular restricted three body problem 
where the inner planet was treated as a test particle. 
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Given the promise of this idea for reconciling planet mi¬ 
gration with the fact that most pairs are not near resonance, 
in this work we extend the theory of eccentricity dependent 
orbital migration and consider some immediate consequences 
of the more complete model. In Section[2] we present simple 
analytic formulae for the stability of the resonant equilibrium. 
These expressions generalize the results previously obtained 
for the circular restricted three body problem. Importantly, 
we show how these criteria depend on the planetary mass ra¬ 
tio and relative eccentricity damping rates between the two 
planets. This allows us to understand for which parameter 
values the instability of the resonance can occur and make 
predictions for which real systems this might have occurred 
for. We test these analytic results numerically in Section R] 
In Section HI we discuss the implications of our work with 
regards to whether or not overstable librations of first order 
resonances can single-handedly account for the observed pe¬ 
riod ratio distribution of Kepler planets. The derivation of the 
formulae presented here is given in the Appendix. 


2 . OVERSTABLE LIBRATIONS OF FIRST ORDER 
MEAN MOTION RESONANCES 

We consider a system of two planets of mass m\ and to 2 
orbiting a star of mass M* with periods near the m:m + 1 
period commensurability. We assume the orbits are nearly 
circular and nearly coplanar. 

In this regime, the Hamiltonian, which governs the conser¬ 
vative dynamics of the planets, is approximately given by 


TT Gtf*rai GM*to 2 Gmim 2 

H =-x 

2ui 2a 2 a 2 


/m+l,27(ares)ei COS [0 - TZ 7i] + 

[/m+1,31 (ctres) l2ct re s]c 2 COS \9 7t7 2 ] 


(1) 


where 6 = (to 4 - 1)A 2 — toAi and a*, A,; and vj % are the 

semimajor axes, eccentricities, mean longitudes, and longi¬ 
tudes of periastron of the two planets. The quantities f rn .+ \ :ij_ 
and fm. 4 -_i si are functions of Laplace coefficients (Murray & 


Dermott 1999 p. 539-556) evaluated at a = <i\ /w> = a res , 
where a res corresponds to exact commensurability. The term 
appearing when to = 1 in the coefficient of the term propor¬ 
tional to e 2 is an indirect term in the disturbing function. 

This Hamiltonian can be reduced to one degree of free- 


dom through a series of canonical transformations (Sessin 
|& Ferraz-Mel lo|p 984[ |Wisdom||1986| |Henrard et al.| 1986| l. 
Since the derivation exists in the literature, we do not repro¬ 
duce it here, though a rough sketch is given in the Appendix. 
After performing the appropriate canonical transformations, 
the Hamiltonian Q. in the region of phase space close to the 
resonance, takes the following form: 


H' = -*($' -T ') 2 


V^COS (j> 


( 2 ) 


where 




1 ( 3.75mA 2/3 

gtJ c 


a 2 sa e\ + e\ — 2eie 2 cos (zui — tu 2 ) 
p = (to + 1) A 2 - \i + ip 


r '~2 


1 /3.75mA 2/3 


c p 




cr z + 


Aa 


tan p = — 


e\ sin vo\ — e 2 sin vj 2 
ei cos tui — e 2 cos vo 2 


(3) 


where Act = a — a res , e p = (mi + to 2 )/M*, p is a gen¬ 
eralized longitude of pericenter, and T' is the “proximity pa¬ 
rameter” which governs how close the system is to resonance. 
All of the above including the Hamiltonian H' are dimension¬ 
less. Additionally, the definitions in Equations ([3]i assume the 
orbits are compact; in this limit, a —> l,m « to + 1 , and 
I/27I ~ /31 ~ 0.8 m. Full expressions, without this “compact 
orbits” approximation, are given in the Appendix. 

As has been pointed out before, the functional form of the 
Hamiltonian in Equation (|2| is identical to that of the circu¬ 
lar restricted three body problem (CR3BP) near a first order 
mean motion resonance, though the parameters have a dif¬ 
ferent meaning. This specific case of the CR3BP has been 
studied at length and functions as a common analytic model 
for res o nances of vari ous type s (se e e.g. Hen rard & Lemaitre 
(|l983j); Murray & Dermott (J1999ji;|Ferraz-Mello[|2007p). It 
is worth noting that in the “compact approximation ’ limit, the 
mass ratio between the planets £ = m\ /m 2 does not appear 
in the Hamiltonian. 

Without dissipation, the proximity parameter is conserved. 
When eccentricities are zero (a = 0), and the orbits are wider 
than the commensurability ( a < a res ),r' is negative, while 
if the orbits are narrow of the commensurability (a > a res ), 
T' is positive. When I - ' < 3/2, the conservative system has 
a single fixed point (x t ), while for T' > 3/2 there are three 
fixed points and a separatrix is present. In this case, two of the 
fixed points are stable (x-\ and x 2 ) and a third (a; 3 ) is unstable. 
To illustrate this, in Figure [I] we show the level curves of the 
Hamiltonian given in Equation Q for two different values of 
the proximity parameter T'. For larger values of T', the unsta¬ 
ble fixed point and the fixed point at the center of resonance 
correspond to approximately equal values of <!>' while the sec¬ 
ond stable fixed point is near zero eccentricity with ( I>' « 0 . 


In the dissipative case, the equations of motion with respect 
to the dimensionless time t' can be written symbolically as 



dH' 

_ 


dt' 

dp 

dt' 

dis 

dip 

dH' 




dt' 

dG_ 

dt 1 


d<&‘ 

dT 

dt' 


dis 


(4) 


In the physical regime of interest, the dissipation acts on 
timescales much longer than the natural timescales of the res¬ 
onance, meaning that the dissipative terms (subscript “dis”) 
are comparatively small in magnitude. We therefore expect 
the system to track on short timescales the level curves of the 
conservative problem with the instantaneous value of T' and 
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Fig. 1.— Level curves of the conservative Hamiltonian for two different 
values of T'. Upper panel: three fixed points as I ' > 3/2, lower panel: a 
single fixed point as T' < 3/2. 


total energy H. 

The extra terms in the equations of motion (|4| are deter¬ 
mined by taking the dissipative evolution of the semimajor 
axes and eccentricities of the planets and converting them into 
the evolution of $' and Trl We parametrize the dissipation 
that acts on the orbit of each individual planet as 

1 de.i 1 

d dt 2 

-w = (- — -—)’ ,5) 

a* at V T e,i T a,i) 

4 At the order of eccentricity we are working, the angle <j>, which depends 
on eccentricities through the generalized pericenter if), does not change due 
to non-conservative effects (see Appendix for proof). 


The above expressions apply for eccentricity decay and mi¬ 
gration towards the star (t Q j) ; and T e ,» are positive; to reverse 
the direction of either the sign of r can be changed). These 
timescales can be estimated analytically (e.g. |Goldreich &| 
|Tremaine||1980[ |Tanaka et al.||2002l |Tanaka & Ward||2004| l 
or from numerical simulations (e.g. |Kley & Nelson|2012 1 or 
taken to be free parameters. Typically r e <A r„ for Type I mi- 
gration of small planets which do not open a gap in the disk. 
In the case of p / 0, there is a coupling between the semi¬ 
major axis evolution and that of the eccentricity evolution, as 
described byjGoldreich & Schlichting|(j2014j). When p = 1, 
the eccentricity damping alone exactly conserves the angu¬ 
lar momentum of each planet, which then decay slowly on 


timescales of r 0 .,, though estimates indicate j> < 1 (Tanaka & 
Ward |20C)4 >. In past studies, a simple exponential damping of 
the eccentricities and semimajor axes was assumed (e.g. |Lee| 


& Peale]2002), without any coupling (p = 0). 

In general, two planets undergoing convergent migration 
will be caught in resonance if the time required to cross the 
resonance width due to migration is sufficiently long com¬ 
pared with the libration period of the resonance and if the 
initial eccentricities of the planets are sufficiently low (e.g. 
Henrard & Lemaitre 1983} Batygin |2015|>. Planets undergo¬ 


ing divergent migration cannot be captured into resonance, 
so throughout this work we focus exclusively on convergent 
migration (r a i > However, even if the inner planet 

originally migrates more quickly than the outer planet, con¬ 
vergent migration can arise after the inner planet reaches the 
inner disk edge, halts and allows the outer planet to catch up. 

2.1. Dissipative dynamics without a-e coupling 

If there is only convergent migration (r fv , —► oc), the dissi¬ 
pative terms are 

dV 


dt 

d4>' 

dt 


dis 


Oo 

T a 


= 0 


( 6 ) 


dis 


where uq i s a po sitive consta nt der ived in the Appendix (see 
Equations ( |A.24| >, Equations ( |A.30[ )) and 


1 

To,2 


1 

T a ,l 


(7) 


For convergent migration, To > 0 and therefore T' grows as 
a —t a res (since er « 0 initially to ensure capture, T' be¬ 
comes less negative). Once the system is caught in resonance, 
a ~ a res and the increase in T' must be compensated with an 
increase in the eccentricities through an increase in a (Equa¬ 
tion 0). Without explicit eccentricity damping, the area en¬ 
closed b y a c ontour of the Hamiltonian is an adiabatic invari¬ 
ant (Henrard;1982). As T' grows, the contour which matches 
the initial area of the trajectory corresponds to larger eccen¬ 
tricities and smaller libration amplitudes. 

When eccentricity damping is included (but p = 0), the 
slow evolution of T' includes an additional damping term, as 
does that of that is, 


dT 

dt 




dis 

dt 


T e 

Co 




(8) 


dis 
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where cq i s a neg ative constant d erived in the Appendix (see 
Equations ( |A.21 1 and Equations ( |A.30| >) which denotes a de¬ 
cay due to eccentricity damping and 



Te,l 7~e,2 


(9) 


centricities as: 


dV 

dt 

d<L>' 

dt 


dis 


dis 


Oq 

T a 


Co + poi 



(13) 


where 


C = rn 1 /m 2 . (10) 

Note first that now the mass ratio between the two planets en¬ 
ters the equations of motion explicitly, while it did not in the 
conservative case or in the conservative case with only migra¬ 
tion. Furthermore, because of the eccentricity damping, there 
is now a value of <f>' = — (ao/co)(r e /r a ) where the change 
in T' is zero. At a fixed migration rate, a shorter timescale 
for eccentricity damping coincides with a smaller equilibrium 
eccentricity since oc a 2 . This equilibrium in a corre¬ 
sponds to an equilibrium of the proximity parameter T' and 
of (I), and hence the fixed point also corresponds to an equi¬ 
librium value of a. For weak dissipation, the equilibrium lies 
close to the stable fixed point in the conservative case (labeled 
x\ in Figure [TJ, assuming a proximity parameter equal to the 
equilibrium value of T' (see Appendix). 

The equilibrium in cr can be turned into an equilibrium 
value for e. \ and e 2 individually by making use of a second 
conserved quantity T 2 , which is zero for circular orbit^j The 
condition that T 2 begins at zero (since the orbits begin circu¬ 
lar due to damping) and remains zero implies that the percen¬ 
ters are always anti-aligned, and that the relationship between 
the equilibrium eccentricities is 


e 2 — C 1 C \/ ^res / A 


( 11 ) 


where 


R = 


|/m+l,27(ctres) I 


[/m-t-1,31 (ctres) 2(5 m ia re s] 


( 12 ) 


The stability of the fixed point for the dissipative prob¬ 
lem can be determined using linear stability analysis. This 
yields three eigenvalues (as we have three dimensions, T', 
and (f >), one of which (ao) is always real and always nega¬ 
tive, along with a complex conjugate pair of the form a± = 
a\ ± ia 2 . a 2 is well approximated as the libration frequency 
of the unperturbed problem evaluated at the new fixed point. 
One finds that a.\ < 0 and so the fixed point is ul tima tely sta¬ 
ble, as has been well established (e.g. !Lee & Peale| <|2002j>). 
The action is no longer an adiabatic invariant, however, and 
the libration amplitude shrinks to zero because the fixed point 
is a stable attractor. 


2.2. Dissipative dynamics with a-e coupling 

|Goldreich & Schlichting|(|2014jl have shown that the stabil¬ 
ity of the fixed point is not guaranteed when the semi-major 
axis evolution is dependent on the eccentricity. The term ap¬ 
pearing when p / 0 changes the dissipative contributions to 
the equations of motion in Equation 0 at lowest order in ec- 


5 Like r', 'I r 2 is conjugate to an angle not appearing in the Hamiltonian, 
and it is therefore conserved. However, unlike T', v l f 2 does not appear as a 
free parameter in the Hamiltonian. See e.g. |Deck et atjf20 1 3j for details. 


where a\ i s a neg ative constant d erived in the Appendix (see 
Equations ( |A.24 1 and Equations ( |A.3Q[ >), not to be confused 
with the semimajor axis of the inner planet, and r a e is defined 


1 _ 1 C tXres 

't~a,e 7"e_ 1 R 'feg 


(14) 


where R is defined in Equation l[T2|) and ( in Equation (fTO}. 


2.2.1. Condition for instability 

Again there is a single fixed point of the dissipative system, 
as in the case where the coupling parameter p = 0. The equi¬ 
librium value for a of the system, determined by the condition 
that dV/dt |di s = 0 and by the relation between and cr, has a 
slight shift compared with the case where p = 0, correspond¬ 
ing to slightly larger or smaller eccentricities depending on if 
r a e is negative or positive, respectively. As before, there are 
three eigenvalues, which we denote as ao and a± = ai±ia 2 . 
As in the case when p = 0, we find that ao is real and always 
negative, so that any motion along the associated eigendirec- 
tion is contracting. 

However, the real part of the complex pair can now be neg¬ 
ative or positive. When 0 | > 0, the fixed point is associated 
with an unstable spiral on the surface spanned by two eigen¬ 
vectors paired with a±. In this case, the lifetime of the system 
in resonance can be finite. The criterion for a\ > 0 is given 
by 


€p ^ €p,crit — 


3 pm 


(1 + C) 2 


623/2 Ta e 


m(C + 1) + Pg^~ 


3/2 


3/2 


(15) 


where e p = (mi + m 2 )/M+, B « 0.8m and r a , r e , and r aj6 

are as defined in Equations 0, 0, and {0- Please note 
that this expression employs the compact approximation ex¬ 
cept for in T a e where we have retained the factors of R and 

ftres¬ 
it is interesting that even in the case where the fixed point 
is associated with an unstable spiral (ai > 0, e p < e PiCr it), 
the evolution still drives the system to the fixed point in reso¬ 
nance before the overall instability drives the oscillation am¬ 
plitude to larger and larger values. Initial capture into reso¬ 
nance requires that the timescale associated with the contract¬ 
ing direction l/|ao| is much smaller than that of 1/ai. The 
evolution of the system near the fixed point is made up of a 
linear combination of the eigenvectors associated with these 
eigenvalues, and as such, on shorter timescales the contract¬ 
ing evolution dominates. On longer timescales, the growing 
evolution takes over, and the system can escape from reso¬ 
nance. This separation of timescales is true especially near 
the critical part of parameter space where the eigenvalue ai 
is changing sign, i.e. where the timescale l/|ai| diverges.We 
will discuss th e tim escales associated with the evolution fur¬ 
ther in Section |4~I1 
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Fig. 2.— Critical curves showing where the fixed point of the 3:2 reso¬ 
nance is associated with an unstable spiral as a function of r ei i/r e> 2 and 
£ = mi /m 2 . The different curves green, black, blue, and purple correspond 
to varying e p , while the red curve shows where r a ^ e = 0. The region above 
the red curve has r a , e < 0 and is stable. Below this curve, the fixed point 
can be unstable, but only if the parameters lie in the lower right region, be¬ 
low the critical curves. The dashed lines show the predictions after making 
the compact approximation (though not for r a , e )- The eccentricity damping 
timescale of the outer planet is fixed at lOOx shorter than the migration time 
Ta, 2, and 7a, i ~ 500r a ,2- 


Regardless of the ultimate stability of the fixed point, the 
equilibrium eccentricity in the resonance is given by 


.2 (1 + C ) 2 

eq 2m(C + 1) + r a ' 


06) 


As expected, if p = 0 and there is no coupling between ec¬ 
centricity damping and semimajor axis evolution, the criterion 
given in Equation © can never be satisfied - the equilibrium 
in resonance is always stable. The equilibrium eccentricity is 
still given by Equation ( fl 6 | . 

A sufficient criterion for stability is e p > e p , cr it- We 
assume that convergent migration leads to resonance cap¬ 
ture and that eccentricity damping leads to an equilibrium 
in resonance. Therefore the parameters r e and r a are posi¬ 
tive. Note that for the equilibrium eccentricity to be real - 
for the equilibrium to even exist - the denominator in Equa¬ 
tion © must be positive, which implies that the factor 
(m((+l)+pT e /Ta : e) 3 ^ 2 appearing in the critical value e PjC7 -it 
for overstable librations is also real and positive. 

This implies that the criterion cannot be satisfied if r a e < 
0, since in that case e p crit is negative. This occurs when 


T e , 2 

T e , 1 


<c 2 


C^res 

~w 


(17) 


In Figure [2] and Figure [3] we show the critical curves gov¬ 
erning the stability of the fixed point for the 3:2 and the 2:1 
resonances, defined by Equation ( fl5j ), on the parameter plane 
°f T e,i/T e ,2 and £ = ) 7 ii/rri 2 . In these plots, we have chosen 
Ta ~ T a 2 ( t 0| i ->■ 00 ) and r e ,2 = r Oj 2 /100. There are no 
other free parameters. 

Each plot shows the following. The red curve denotes 
T a, e = 0 , and the area above the red curve corresponds to 
a region of parameter space where the fixed point is stable 
since r a , e , < 0. Below the red curve, the fixed point may 


FIG. 3.— Critical curves showing where the fixed point of the 2:1 resonance 
is associated with an unstable spiral as a function of T e i/r e 2 and £ = 
m\/m 2 - Refer to Figure[2]for details. 


be unstable, but only if e p < e PiCr it. The green, black, 
blue, and purple curves correspond to e p = f p . crl t for e p = 
( 10 - 3 , 10 ~ 4 , 10 — 5 , 10 ~ 6 ), respectively, and the area to the 
lower right of these curves is where the instability occurs. 

In order to simplify the expressions for the critical values 
of e.g. e p to that given in Equation ( p~5] >, we made the “com¬ 
pact approximation” that a = ai/a 2 — > 1 as discussed above 
(note again that we do not make that approximation for r a , e ). 
This approximation is poorest for the 2:1 resonance, both be¬ 
cause o: res is further from unity but also because of the indirect 
contribution to the coefficient of the e 2 term in the disturbing 
function, so that R = 2.78. However for closer resonances 
the approximation is very good. 

The difference between the dashed (approximate) and solid 
(exact) e p = e PiCr it curves in the Figures demonstrates this. 
Note that the disagreement would be stronger for the 2:1 res¬ 
onance if we had used the compact approximation for r o e , 
which is why we keep the full expression (it is easy to include 
it here and retain a simple expression for e PtCr iu the same is 
not true if we had never made the approximation at all). Re¬ 
gardless, the exact formulae is reasonably well approximated 
by the estimate, and both show the basic result that systems 
with a more massive planet interior are more stable against 
overstable librations. This is especially striking when we con¬ 
sider the two cases of an inner and an outer test particle below. 

2.2.2. Limiting case of CR3BP 

Given that these results were derived within the framework 
of the elliptic planetary three body problem, they should re¬ 
duce to the results obtained previously in the limit of the cir¬ 
cular restricted three body problem. We first consider the case 
of a test particle moving outwards towards a massive planet. 
In this case, —>• 0 , r a , 2 —> 00, r e , 2 —> 00, and r 0jl —» — t 0|1 
(to account for outward migration). Then r e = r a , e = r e ,i 
and T a = —T a . i- We also assume the outer planet has a circu¬ 
lar orbit, so a = e\. In this case, when p = 1 the equilibrium 
eccentricity is given by 


2(m + 1) r a ,i 3(m + 1) r„,i 

(18) 
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where r n _\ = 2/3r a ,i . This agrees with Equation (24) of Gol- 


dreich & Schlichting (2014 1 . The criterion for overstabihty is 

3/2 


£2 < 


m 


_ ( 

B v / 3(m+ 1)3/2 \ Tn l 


(19) 


which agrees with Equation (30) of Goldreich & Schlichting 
(20140 

We now turn to the opposite case of a test particle mov¬ 
ing inwards towards a massive planet. In this case, f —► 
oo, r a> i -A oo and r e> i -A oo. Then r e = r e , 2 /C, t q = r 0)2 , 
and r o e = —Te^/C, 2 - The fact that r o e < 0 in this case im¬ 
mediately implies that the resonance is stable. Indeed, after 
carefully taking the limit as C, > oo, the criterion for over 
stability is 


3 pm 1 

623/2 ( m — p) 3/ 2 



( 20 ) 


which can never be satisfied. When m < p, this expression ei¬ 
ther diverges or becomes imaginary. However, in these cases, 
the equilibrium itself does not exist (see Appendix). Note also 
that the divergence when m = p = 1 does not occur within 
the full expression (i.e. without taking the compact limit). 

Why is there a difference between the two limiting cases? 
We do not yet have a good physical intuition for this. How¬ 
ever, since the Hamiltonian itself is approximately indepen¬ 
dent of the mass ratio between the two planets, any depen¬ 
dence on the mass ratio must come from the dissipative terms. 
In the case where the inner planet becomes a test particle, the 
eccentricity damping leads to an inward migration of the test 
particle proportional to e\, in opposition to the overall out¬ 
ward migration towards the massive planet. When the outer 
planet is the test particle, the eccentricity dependent migra¬ 
tion acts coherently with the direct semimajor axis damping 
to move the test particle towards the inner planet. The small 
contribution of p in one case apparently compromises the sta¬ 
bility of the resonance, while in the other it stabilizes it fur¬ 
ther. 


2.2.3. Condition for escape from resonance 

If the fixed point is an unstable spiral, the libration ampli¬ 
tude about the fixed point will grow in time. For a range of 
parameter values, these oscillations will saturate at a stable 
limit cycle enclosing the unstable fixed point. In this case, the 
system remains trapped in resonance but with a nonzero li¬ 
bration amplitude. A criterion for escape from the resonance 
would be such that the fixed point is unstable and there is no 
possibility of saturating at a stable limit cycle. 

We can understand this qualitatively as follows. The true 
fixed point of the dissipative problem lies near the fixed point 
of the conservative problem at the center of the resonance re¬ 
gion (labeled xi in Figure [lj. There exists a limit cycle be¬ 
cause the motion is being driven towards the fixed point in one 
eigendirection and away from the fixed point as an unstable 
spiral in the other two eigendirections, and there is a balance 
of these opposing actions at some point. However, when there 
are three fixed points of the conservative problem (T' > 3/2), 

6 In their formulation, they use x„ rather than r a . Additionally, though 
we have used the same symbol for the coupling parameter p, our case with 
p = 1 corresponds to their case with p = 3. Finally, we are using a capital 
B to represent 0.8m, which they use a lowercase ft for, because we use a 
lowercase 0 for a different meaning in the derivation in the Appendix. 


and the fixed point of the full dissipative problem is an unsta¬ 
ble spiral, the libration amplitude will grow until the trajectory 
enters a region of attraction in the (v / 24>' cos <f>, v / 24>' sin </>) 
plane near the conservative fixed point cc 2 at a ss 0, without 
reaching a stable limit cycle. 

This region of attraction corresponds to what would have 
been the inner circulation region (corresponding to oscilla¬ 
tions about x 2 ). The point X 2 is not a fixed point of the 


full problem, but the motion in the (v / 2$' cos <j>, v/^Nsin (j>) 
plane near this region appears as a spiral towards X 2 , because 
the dissipative evolution approximately follows contours of 
the conservative problem on short timescales. For large values 
of T', the stable fixed point (x 2 ) of the conservative problem 
corresponds to nearly zero eccentricity ct, and so in this at¬ 
tractive region the eccentricities of the planets damp to nearly 
zero eccentricity. T' continues to grow, since this is not a 
fixed point of the dissipative problem, and this brings the ec¬ 
centricities closer to zero and brings the pair narrow of the 
resonance (since for o ~ 0 a positive T' implies a > a res , 
see Equation 0)- This evolution is illustrated in Section 13.1 


where we show the numerically determined evolution ot 


and <f> on the (v / 2$' cos <j), y/2^' sin <j>) plane along with ap¬ 
propriate contours of the conservative Hamiltonian. 

Applying this criterion, we find that the system avoids being 
trapped in a limit cycle if e p < e PjCr n and T' > 3/2 or 


e P ;$ 


3 TO 2 


(1 + C) 3 


166\/2 / \ 3/2 

( m(C + 1 )+PTff 


3/2 


_ rn r a , e (1 + C) , 

^LC — q €p,crit 


= £LC 


( 21 ) 


If 1) £p,crit < c P or 2) e LC < e P < Ep, cr it the system is 
stuck in resonance, either at the (stable) fixed point in the for¬ 
mer case or in a limit cycle about the (unstable) fixed point in 
the latter. Note that the critical value epc is nonzero even 
if the coupling parameter p = 0. This doesn't mean that 
the system can escape from resonance if e p < e rc even if 
p = 0, but that this limit cycle criterion is meaningless un- 
less e p < e v ,crit in the first place. Finally, as pointed out in 


Goldreich & Schlichting (2Q14|>, the limit cycle is only a fac¬ 


tor for resonances satisfying epc < ( pcrlt . For an inner test 
particle with p = 1, this corresponds to m < 8. Because of 
this, we focus on Equation ( fl5] l as a criterion for instability, 
though formally one requires e p < elc < fp.c.ni for escape 
from resonance. 


3. NUMERICAF RESULTS 

Here we test how well our simple analytic criterion applies 
to a “real” system using direct numerical integration of the 
full gravitational equations of motion with the appropriate 
dissipative terms put in. We integrate the standard gravita¬ 
tional equations of motion using a Bulirsch-Stoer integration 
scheme. The migration terms are added directly to the equa¬ 
tions of motion following the prescription in the Appendix of 
jLee & Peale ( |2002| >. This requires applying the chain rule 
to determine how changes in a and e, defined in Equation ([5]), 
translate into changes in the cartesian positions and velocities. 
As the planets migrate towards the host star, their orbital pe¬ 
riods decrease, and an adequate time step for the initial orbits 
may be too large for the orbits at a later time in the integra¬ 
tion. To alleviate this issue and keep the (fixed) dissipation 
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timescales slow compared to the orbital periods, we rescale 
the semimajor axes of the planets at each time step so that the 
semimajor axis of the inner planet is fixed. Then our default 
time step is always short compared to the orbital periods, and 
the fixed migration rates are typically long compared to the 
the relevant libration timescales. 

3.1. Illustration of instability effect 

We begin by showing the explicit evolution of 0 , and 
T' in the three cases of permanent capture with no limit cy¬ 
cle, permanent capture with a limit cycle, and escape from 
resonance to better illustrate the above discussion of Section 

12331 

In Figure H] Figure [5] and Figure [ 6 ] we show the evolution 
of three different systems. We have set T ei 1 = r e , 2 , C = 0 . 6 , 
f e ,2 = 10 4 Pi, and T a .i r a 2 - In case 1 e p = 8 x 10 -4 and 
T a ,2 ~ 10 6 Pl, in case 2 e p = 8 x 10 -5 and r a> 2 ~ 2 x 10 6 Pi, 
and in case 3 e p = 8 x 10 5 and r a 2 ~ 10 6 Pi. In case 1 and 
case 3, e p , cr it = 3 x 10 -4 and £lc = 10 -4 , while in case 
2 tp,crit. = 10 - 4 and clc = 3.5 x 10~ 5 . To calculate these 
values, we used Equations ( fT5j ) and ( |2Tj ). 

In all cases, the planetary system is captured into the 2:1 
mean motion resonanc e and i nitially driven to the fixed point, 
as discussed in Section |2.2.1| We show in the upper panel of 
each of these figures the contours of the conservative Hamil¬ 
tonian on the (x/2<3>' cos </>, sin ff) plane (in red) corre¬ 
sponding to approximate equilibrium value of F'. Overplot¬ 
ted on these contours is the behavior of the system variables 
undergoing full dissipative evolution. In all three cases, the 
system begins at zero eccentricity (the origin). The planets 
are captured into the mean motion resonance and this initially 
leads to an increase in eccentricities (radial distance from the 
origin) until the equilibrium is reached (this stage of the evo¬ 
lution is shown with cyan points). 

First, in Figure [4] we show case 1, where we have chosen 
e p > £p,crit so that the fixed point is stable (a\ < 0). The sys¬ 
tem remains at the fixed point. In the bottom plot, we show the 
evolution of F', <!>', and the fractional deviation in the period 
ratio from the exact commensurability. It is difficult to see by 
eye, but the amplitude of oscillation of $' is decreasing as we 
would expect. 

In Figure B] we show case 2, where we have chosen clc < 
fp < c 7J ,. r)t 3’hc limit cycle behavior is possible because the 
equilibrium T' is less than 3/2, and so there is only one fixed 
point of the conservative problem. In the upper panel we now 
show in dark blue the stage of the evolution where the am¬ 
plitude of oscillations grows about the fixed point (instabil¬ 
ity). In grey the limit cycle is shown. Note that it encloses 
the origin, and therefore the resonant angles are circulating in 
this configuration. The amplitude of the limit cycle is chang¬ 
ing because the proximity parameter l 7 is oscillating as well, 
which changes the level curves of the conservative Hamilto¬ 
nian. It is unclear what exactly causes the behavior at a time of 
0 . 4 T a , when the amplitude of oscillations suddenly decreases 
and then begins to increase again. We note that a precise an¬ 
alytic description of the limit cycle is quite complicated; the 
limit cycle criterion we use is a heuristic one. 

Finally, in Figure[ 6 ] we show how a system can escape from 
resonance. In this case, e p < e^c < £p,crit because T' is 
greater than 3/2 in the equilibrium configuration. The ampli¬ 
tude of oscillation grows about the fixed point (shown again 
in blue) until it crosses into the basin of attraction dominated 
by the second stable fixed point of the conservative problem. 



(20’) cos[(|)] 



0.0 0.1 0.2 0.3 0.4 0.5 


Fig. 4.— Permanent capture into the 2:1 resonance in the case where the 
fixed point is stable. Upper panel: contours of the conservative Hamiltonian 
at the equilibrium value of H Ri —0.52 (red) and the actual dissipative evo¬ 
lution of the system showing capture into resonance (cyan). Lower panel: 
time evolution of E (red), E (black), and the fractional deviation of the pe¬ 
riod ratio from 2.0 in percent (blue). One can see a small damping of the 
oscillation amplitude of <E>'. 


This stage of the evolution is shown in black. At this point, 
begins to decrease while I’ 7 continues to increase (lower 
panel). As V increases, the second stable fixed point of the 
conservative problem moves closer to the origin and the sys¬ 
tem eccentricity decreases to zero. This then causes the period 
ratio to shift to be narrow of the resonance since F' relates the 
eccentricity er and the period ratio. 

3.2. Long-term evolution 

In Figure [7] we show the evolution of the period ratio and 
eccentricity a of two “restricted-like” systems. On the left 
side, the two plots show the results when the inner planet 
is effectively a test particle with a mass 10~ 8 the mass of 
the star, while the outer planet is HE 5 the mass of the star. 
We choose the migration rate and eccentricity damping of the 
outer planet to be very long compared with all other physical 
timescales. The migration timescale of the inner “test parti¬ 
cle” outwards is approximately 10 6 times the orbital period 
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Fig. 5.— Permanent capture into the 2:1 resonance in the case where the 
fixed point is unstable. Upper panel: contours of the conservative Hamilto¬ 
nian at the equilibrium value of T' « 0.49 (red) and the actual dissipative 
evolution of the system showing capture into resonance (cyan) and subse¬ 
quent growth of the libration amplitude (blue) to a stable limit cycle (grey). 
We do not show the entire evolution from initial onset of instability to the 
limit cycle. Lower panel: time evolution of T' (red), d?' (black), and the frac¬ 
tional deviation of the period ratio from 2.0 in percent (blue). The dashed 
lines and the regions in between reflect those used for the “capture”, “onset 
of instability”, and “escape” evolution in the upper panel. The sudden change 
in amplitude of oscillations at ~ 0.4r a is not captured by our simple analytic 
formulation. 


of the test particle and the eccentricity damping timescale is 
chosen to be lOOx smaller. Both bodies begin with circular 
orbits outside the 2:1 resonance. In this case, the criterion for 
whether or not the instability can set in at the 2:1 resonance 
is given by TO 2 /M* < 5 x 10“ 4 , which is easily satisfied (re¬ 
gardless of whether we use the compact approximation or the 
exact form). Indeed, we see this single trajectory undergoing 
capture and then subsequent escape via overstable librations 
for many resonances. The dotted colored lines show the pre¬ 
dictions for the equilibrium eccentricity of the test particle 
for each resonance (without making the approximation that 
R « 1 or to ~ to + 1). 

In the right set of panels, we show the results of integra¬ 
tions of a series of systems where the outer planet is treated 


Fig. 6.— Temporary capture into the 2:1 resonance in the case where the 
fixed point is unstable. Upper panel: contours of the conservative Hamilto¬ 
nian at the equilibrium value of T' rj 1.86 (red) and the actual dissipative 
evolution of the system showing capture into resonance (cyan), subsequent 
growth of the libration amplitude (blue), and finally escape from the reso¬ 
nance by damping to second stable fixed point of the conservative problem 
(black). See text for an explanation. Lower panel: time evolution of T' (red), 
(black), and the fractional deviation of the period ratio from 2.0 in percent 
(blue). The dashed lines and the regions in between reflect those used for the 
"capture”, “instability”, and “escape” evolution in the upper panel. 


as a test particle. In this case, we reversed everything exactly 
compared with the left set of panels. The outer “test particle” 
has a mass of 10 -8 the mass of the star, while the inner planet 
is 10~ 5 the mass of the star. The migration and damping rate 
for the inner planet is very long compared with all other phys¬ 
ical timescales. The migration timescale of the outer test par¬ 
ticle is again approximately 10 6 times the orbital period of the 
test particle and the eccentricity damping timescale is chosen 
to be lOOx smaller. 

In the upper right plot, the red curve shows the result when 
both bodies begin with circular orbits outside the 2:1 reso¬ 
nance, the other colors show the evolution when began just 
outside of other first order mean motion resonances as labeled. 
For each, the system is captured into resonance and does not 
escape. That is, the libration amplitude does not grow with 
time, suggesting that these are in fact stable configurations. If 
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Fig. 7.— The asymmetry between an outer test particle and an inner test 
particle. On the left we show the evolution of a single trajectory in period 
ratio (top) and eccentricity (bottom) of a test particle migrating outwards to¬ 
wards a massive planet on a fixed circular orbit. On the right we show the 
evolution of four separate orbits, began just outside the 2:1 (red), 3:2 (blue), 
4:3 (green), and 5:4 (purple) resonances, for an outer test particle migrat¬ 
ing inwards towards a massive planet on a fixed circular orbit. The dashed 
colored lines in the bottom two panels show the estimated equilibrium eccen¬ 
tricity of each resonance from the exact expression derived in the Appendix. 
See text for further details and discussion. 


3:2 MMR 



Fig. 8.— The critical curve (black) showing where overstable librations 
of the 3:2 resonance can occur for a system of two planets of total mass 
e p = 1.31 X 10~ 5 . The evolution of the three systems shown with constant 
(j (filled circles) is shown in Figure[9] while the evolution of the fou r sy stems 
shown with constant T e ,i/r e ,2 (open circles) is shown in Figure MO The 
color of th e points corresponds to the color of the trajectories in Figure 9|and 
FigureflO] 


we took the results for the inner test particle case and applied 
them here, we would expect these systems to easily be unsta¬ 
ble as well since my/M* < 5 x 10~ 4 . Moreover, if the outer 
test particle case was analogous to the inner test particle case, 
the timescale for escape would be much shorter than our inte¬ 
gration time, as it was in the case when the inner planet was a 
test particle. This increases our confidence that we have inte¬ 
grated these systems long enough to show that they are stable 
(i.e. that the instability time is not significantly longer than 
the integration time, giving the illusion of stability), and that 
there truly is a difference between whether the inner planet is 
a test particle or the outer planet is a test particle. 

We now turn to some tests of the criterion for compara¬ 
ble mass planets. We set T a ^/Pi ~ 10 6 , 744 = 5 r aj 2 , 
and r e 2 = t Oj 2 / 100 . We then varied r e>1 and the ratio of 
£ = ini/rri 2 ■ In Figure [I] we show on a panel of r ei i/r ei 2 
and mi/m 2 where overstable librations can occur for these 
parameters. The points on this plot correspond to systems we 
studied migrating convergently beginning from just outside 
the 3:2 resonance. We show the evolution of the three sys¬ 
tems with a fixed mass ratio £ = 10 -05 (marked with filled 
circles in Figure [ 8 | in Figure [9] As predicted, the two sys¬ 
tems with larger values of T ei i/T e ,2 (those in red and blue) 
undergo instability and escape from resonance on a timescale 
short compared to the integration time, while the system with 
the smallest value of r ei i/r e ,2 (in black) remains in resonance. 
Note that all of these systems have a value of r a e > 0. 

In Figu re[l0| we show the evolution of four systems marked 
in Figure [ 8 ] with open circles with r e ,i = r e 2 but varying 
mass ratio ln this case, systems with £ > 1 should be stable 
against overstable librations because the quantity r aj£ . is neg¬ 
ative. This is indeed what is observed. If £ < 1, the orbits 
may be unstable, but only if the total mass of the planets is 
low enough. For the orbits shown in purple and black, that re¬ 
quirement is satisfied, as shown in Figure [ 8 ] (the open circles 
in purple and black lie below the solid black curve), and these 
orbits do escape from resonance. 


1.6 
.2 1.5 
m j 4 
DC 1 ' 
T3 1.3 

O J _ 
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Fig. 9.— The evolution of three systems with e p = 1.31 X 10 -5 , 
C = 10 -0 ' 5 , t 0 , 2 /Pi ~ 10 6 , t 0 ,i = 5r a , 2 , and T e ^ 2 = T a , 2 /100. 

These trajectories, and their color, correspond to the three filled dots shown 
at constant C in Figure[8] At this mass ratio, as predicted, only systems with 
T e i > 10 — °' 8 r e ,2 = 0.15r ej 2 undergo overstable librations. 
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4. DISCUSSION 
4.1. Timescales for escape 

The timescale on which the pair of planets escapes from res¬ 
onance after reaching an equilibrium eccentricity larger than 
the critical value is an important quantity because when com¬ 
pared with the migration time r a it will determine what frac¬ 
tion of the time any given pair is found in a resonance vs. 
migrating between them. The evolution shown in Figure [7] 
Figure [9jand Figure p~0] indie ate that the systems which do es¬ 
cape from resonance typically spend the majority of their time 
trapped in resonances undergoing overstable librations rather 
than in between them. 

The real (positive) part of the eigenvalue au gives the rate 
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Fig. 10.— The evolution of four systems with e v = 1.31 X 10 5 , 
T a ,l/P\ ~ 10 6 , Ta,i = 5 t Bi2 , T e ,2 = T a , 2 /100, and r e ,i = r e , 2 . These 
trajectories, and their coloncorrespond to the four open circles shown at con¬ 
stant r e ,i/r e ,2 in Figure [8] At this ratio, as predicted, only systems with 
mi < m 2 undergo overstable librations. (note that the red curve lies below 
the blue curve in terms of period ratio evolution in the upper panel; the two 
are distinct in the lower plot showing equilibrium eccentricity, however). 


at which orbits spiral away from the fixed point when the or¬ 
bits are near the fixed point (as it is a local stability analysis). 
Since the real part of the eigenvalue a \ passes from negative 
to positive on the critical curve e c = e pj . r ,t, the instability 
time, given by the inverse of the eigenvalue, is very long near 
the critical curve. The instability time is not the same order 
of magnitude as r e because the eigenvalue is proportional to 
1 /t,. and the difference of comparable quantities (unless the 
system satisfies e p cr it) or vice versa). On the other 

hand, the timescale associated with the negative eigenvalue 
ao, which is also proportional to l/r e , is much shorter than 
timescale associated with the positive eigenvalue a± since it 
does not have this dependence on a difference between like 
quantities. In fact, temporary capture necessitates that the in¬ 
stability timescale be much longer than ~ 1 /1 cko | - If the two 
timescales were comparable, the system would not even be 
captured into resonance in the first place. This implies that 
escape from resonance must occur on timescales significantly 
longer than r e in general. 

To illustrate these points, we show in Figure m the 
timescales r e , a)” 1 , and |ao| _1 , in units of r a , for the 3:2 
resonance for a range of r e i/r ej 2 - The inner orbital period 
is ss 10 2 days. All parameter values correspond to those of 
Figure [I] except ( is fixed at 10 -0 ' 5 . One can see that a) 1 
(shown only when positive) diverges when the critical curve is 
reached at Log[r ei i/r ej 2 ] ~ — 0.8 and at Log[r ej i/r ej 2 ] ~ 0.9 
as in FigureJS] Moreover, across the entire range, r e is compa¬ 
rable to | ao r 1 and at least an order of magnitude smaller than 
a)” 1 . The numerical simulations show that the pairs spend 
most of their time in resonance, which suggests that approxi¬ 
mately 5-10 instability times are required for escape (the for¬ 
mer implies roughly T escape ~ r a and Figure[TT]indicates that 
r Q ~ 10a)” 1 , and therefore T escape ~ 10a)” 1 


3:2 MMR 



Fig. 11.— Relevant timescales to the problem of dissipative evolution near 
the 3:2 mean motion resonance. Here we show r e and the two timescales 
associated with the real eigenvalues determined by a linear stability analysis 
near the fixed point in units of r a . The eigen-timescale l/cx\ is only shown 
when cx.\ > 0 and the fixed point is ultimately unstable. The relevant param¬ 
eters used are given in the text. 


served period ratio distribution because 1) the typical masses 
of the Kepler planets are small enough, given the estimated ra¬ 
tio of migration and eccentricity damping rates expected from 
Type I migration, that the first order resonances are unstable 
and any capture was only temporary, and 2) the time spent in 
resonance was small compared to the time spent in between 
resonances, such that when planet-disk interactions stopped, 
most pairs were left in between resonances. At the same time, 
this theory is consistent with giant planets being found in res¬ 
onance ( jWright et al.|20lf) , since in that case the total mass 
of the planets is too large compared to the critical value for 
instability. 

In the more general case of two massive planets, we have 
shown that overstable librations can only occur in specific 
cases. First, the resonance is stable if the quantity T a e is neg¬ 
ative, regardless of the total mass of the planets or the ratio of 
eccentricity damping to semimajor axis damping, that is, the 
resonance is stable if 

r e,2 rn 'i , jj2 

^ 2 ^Tes/rt (22) 

T e ,l TO 2 


where R and a res are of order unity. If we employ the “com¬ 
pact approximation” and the scaling that t p oc where p 
stands for planet and ra_ damping timescale due to interac¬ 


tion with the disk (e.g. Goldreich & Tremaine (1980); Tanaka 


et al.|(|2002|i), we find that resonances are nominally stable to 
overstable librations if 


4.2. Planet pairs discovered with Kepler 

According to jGoldreich & Schlichting ( j2014) >, overstable 
librations of first order resonances could account for the ob¬ 


777.2 < 7771 (23) 

Note that for the 2:1 resonance there is considerably more 
leeway in this mass criterion, as now the pre-factor on the 
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right-hand side of Equation (|22| evaluates to ss 1/12. For 
comparison, at the 3:2 resonance it is 1/0.85. This implies 
that if r p (X —, the 2:1 resonance is stable for 


m 2 < mi/12, 


(24) 


regardless of the total mass of the planets. 

If the system does fail Equation ( 22|i, we find a similar result 


Goldreich & Schlichting 


That is, the system 


to that of 

can escape from resonance via overstable librations only if 
the total mass of the planets is small compared to a critical 
value. The critical value depends on the ratio of r a /r e , ( and 
the relative eccentricity damping times. Therefore, a pair con¬ 
taining a giant planet (regardless of planet mass ratio) will be 
more likely to be caught in resonance permanently given the 
high value of e p , while low mass planet pairs like those found 
by Kepler will be liable to escape from resonance. What this 
means is that escape from resonance should have occurred for 
a subset of the Kepler systems. In particular, the lack of pairs 
caught in the 2:1 mean motion resonance can be explained 
through overstable librations assuming the total mass of the 
planets is small enough and that m 2 < mi/12. However, 
unless the inner planet is less massive than the outer planet, 
capture into all other first order resonances is predicted to be 
permanent. 

More importantly, even for pairs which undergo overstable 
librations, the cumulative time that the pair spends in any res¬ 
onance is much longer than that spent in between resonances. 
Therefore the average period ratio distribution of a popula¬ 
tion of pairs undergoing convergent migration will be peaked 
near resonances, since all pairs will either 1) be too massive 
to undergo instability or not satisfy the mass ratio £ require¬ 
ment (mi < m 2 ), and therefore be captured permanently into 
resonance or 2) satisfy the total mass and mass ratio require¬ 
ment for instability but spend most of their time in a resonance 
regardless. Overstable librations alone cannot explain the pe¬ 
riod ratio distribution. 

One possible resolution to this issue could be orbital insta¬ 
bility. If any pair in a system satisfies the requirements for 
overstability, that pair will undergo successive repetitions of 
the following sequence: capture into resonance, followed by 
overstable librations and escape from resonance. This would 
lead a gradual compactification of the orbits - and could ex¬ 
plain the presence of some very close pairs of planets, like 
Kepler-36 (Carter et al. 2012| >. At some point, if this process 
continues, the pair will be close enough that when they es¬ 
cape the current resonance, at the equilibrium eccentricity of 
the resonance, they are in the chaotic region corresponding 
to overlap of resonances (e.g. [Wisdom] (1980); Deck et ah 
(|2013|). If the timescale to develop crossing orbits in this 
cEaotic region is short enough, these planets will scatter or 
collide with each other. If they do not collide, and merely 
scatter, they will begin to undergo the entire process again, as 
long as the disk is present. 

When the disk finally dissipates, the planetary system will 
be left in a precarious position, with compact pairs of plan¬ 
ets in resonance with nonzero eccentricities of several percent 
(corresponding to the equilibria at the center of resonance). 
Pairs that were undergoing overstable librations will have 
nonzero libration amplitudes. We hypothesize that these con¬ 
figurations will typically have lifetimes short compared with 
the ages of the systems, and that many of them will undergo 
orbital instability, leading to wider distribution of period ra¬ 
tios and pairs not near resonance, in correspondence with the 


period ratio distribution. This would also explain why there is 
a decrease in the number of very compact pairs as well, since 
those orbits would be most unstable (see also |Pu & Wu|2015j ). 

4.3. Effects of external sources of precession 

The analysis we have undertaken applies when the mo¬ 
tion of the planets is governed by the Hamiltonian given in 
Equation (JTJ and the dissipation prescription of Equations ([5]). 
This formulation neglects apsidal regression of the pericenters 
induced by the gravitational potential of the disk (Heppen-| 
heimer 1980. Tamay oet al.|2015 1 and precession of the peri¬ 
centers directly caused by the density wa ves excited by the 
planet in the disk ( Tanaka & Ward||2004) . If these “forced” 
precession rates differ between the two pl anets, the conserva- 
tive Hamiltonian is no longer integrable (El Moutamid et al. 
2014). Instead, both resonant angles fff = 0 — w 1 and 
Uf = 9 — vj -2 are independent. The resonant centers, where 9 1 
and 62 have zero time derivative, are separated by an amount 
of w -2 — vj\. This splitting changes the period ratio which 
corresponds to the resonance as P2/P1 = (m + l)/m — 
ztiP 2 /{2'Km). 

Although a study of the effect of these external precession 
rates is beyond the scope of this paper, we can provide some 
general hypotheses as to how they might affect our results. 
First, in the limit of a = < 21/02 —> 1, the precession rates 
caused by the gravitational potential of the disk - which are 
independent of the planet masses and but depend on the semi¬ 
major axis of the planets - will be approximately equal. If 
these dominate the external precession rates, there will be no 
splitting between the resonant angles 9 \ and 62, and the inte¬ 
grable analysis should still apply. 

Second, if the precession rates cause a large enough split¬ 
ting between the two resonances and they do not perturb each 
other too strongly, they can be treated individually. In this 
case, the motion would still be integrable near each resonance. 
If, for example, the system was near the resonance associated 
with 9 \, we could average the Hamiltonian over the angle 6 * 2 , 
and we would be left with something akin to the circular re¬ 
stricted three body problem for an inner test particle. In the 
opposite case, we would be left with a problem akin to the 
circular restricted three body problem for an outer test parti¬ 
cle. As demonstrated by|Goldreich & Schlichting|(|2014], the 
CR3BP with an inner test particle is susceptible to overstable 
oscillations, and as shown here, the opposite is not. 

As the system evolves, the period ratio will approach the 
value (m + l)/m from a larger value. The system will first 
encounter one of the resonances (which it encounters first de¬ 
pends on the precession rates). If the migration rate is slow 
enough and the eccentricities low enough, the system will be 
captured. If this first resonance is the O2 one, the system will 
be trapped permanently. However, if the system first encoun¬ 
ters the 61 resonance, the capture may be temporary if the sys¬ 
tem satisfies the criterion on e p . The system could then escape 
the resonance - but it would then encounter the 62 resonance. 
Any capture again would likely be permanent. 

If the splitting between the two resonances is not large 
enough to treat them individua lly, the motion may be ch aotic, 
due to overlap of resonances (El Moutamid et al. |2014] >. In 
this instance the probability of capture and of escape via over¬ 
stable librations might be best assessed through numerical ex¬ 
periments. 

5. CONCLUSION 
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In order to understand the observed sample of exoplanets, 
we need a better understanding of how various physical pro¬ 
cesses shape planetar y sys tems after formation. We have built 
upon the work of |Goldreich & Schlichting] < )2Q 14) > to account 
for the effects of eccentricity dependent serminajor axis evo¬ 
lution for a system of two massive planets with mildly eccen¬ 
tric orbits undergoing convergent migration near a first order 
mean motion resonance. 

We have shown that in the full problem resonances are in 
general more stable against overstable oscillations compared 


with the restricted case considered by Goldreich & Schlicht- 
(f2014jl. 


ing 
than the 


In particular, for first order resonances other 


resonance, if r ei 2 /r ei i < resonance 

capture is permanent (sans other disruptive effects like tur¬ 
bulence), regardless of the total mass of the planets. If the 
timescale for damping is inversely proportional to the mass 
of the planet, then this criterion states that the resonances are 
stable if m 2 < m\. For the 2:1 resonance, however, the same 
is true only if the inner planet is significantly (~ 12x) more 
massive than the outer. 

Now, if the relative eccentricity damping rates and the mass 
ratio between the planets do satisfy the criterion T ej 2 /T e 1 > 
m i/ m 2 ’ we a resu lt ver y similar to that found in the re¬ 
stricted case: resonances are unstable if the total mass of the 
planets is smaller than a critical value dependent on the gen¬ 
eralized r a and r e defined above. This is consistent with the 
fact that pairs of gas giants have been found in resonance (e.g. 


| Wright et al.fcOll) . 

We have also demonstrated that in the case where the res¬ 
onances are unstable the timescale to escape from resonance 
is not simply proportional to r e <C t„. Instead, the instability 
timescale is much longer than r e , and furthermore ~5-10 in¬ 
stability times are required for escape. Because of this, a pair 
of planets undergoing overstable oscillations spends the ma¬ 
jority of time in a resonance rather than in migrating between 
them. 


It is important to determine if the lack of period commen- 
surabilities within the Kepler data is evidence that convergent 
migration did not occur. This paper has attempted to address 
whether the mechanism proposed by (Goldreich & Schlichting| 
(2014) can explain the lack of near resonant pairs within the 
context of orbital migration. Testing predictions about which 
Kepler pairs could have escaped from resonance based on 
their mass ratio is difficult since the Kepler planets do not typ¬ 
ically have measured masses, and the radius ratio may not be a 
reliable proxy ( |Weiss & Marcy|2014[ >. However, based on the 
timescale argument, even if the majority of pairs did satisfy 
the criterion for overstability, the lack of near resonant pairs 
within the Kepler period ratio distribution is not due to pairs 
simply spending more time between resonances than trapped 
in them, so that when the gas disk dissipated more pairs were 
left between resonances than near them. From this we con¬ 
clude that overstable librations, which should have occurred 
for some pairs, cannot single-handedly account for the lack of 
pairs near resonance. 

We have argued that regardless of whether planets form in- 
situ or at larger orbital distances, they likely interacted with 
a gaseous disk, the presence of which may be required to 
explain the volatile rich nature of the exoplanets larger than 
1 .6i?0. In this case, then, a loss of orbital energy and some 
exchange of angular momentum with the disk is unavoidable, 
and therefore migration must be brought into agreement with 
the lack of near resonant pairs. An interesting idea left for 
further investigation is that the effects of higher multiplicity 
systems and of true orbital instability could account for the 
period ratio distribution in combination with overstable libra¬ 
tions. 


K.M.D. would like to acknowledge support from the Joint 
Center for Planetary Astronomy at Caltech and also to thank 
Peter Goldreich for prompting the discussion of external 
sources of precession. 
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APPENDIX 


A-l. Setup of integrable problem 

We consider a system of two planets of mass mi and m 2 orbiting a star of mass M* with periods near the m:m + 1 first 
order mean motion resonance. We assume the orbits are nearly circular and nearly coplanar, and at this stage we ignore external 
sources of pericenter precession which arise from the gravitational potential of the protoplanetary disk as well as from direct 
interactions with the disk. In this configuration, the Hamiltonian can be reduced to a one-degree of freedom system with a single 
free parameter (Sessin & Ferraz-Mellojl984jl. This is because, after expanding the Hamiltonian about the resonance center, there 
is a series of canonical transformations which reduce the number of degrees of freedom from 4 to 1 dWisdom|1986||Henrard et al.l 
|T986}. 

Before presenting that simplified Hamiltonian, we motivate it as follows. As given in Equation (jJJ), the Hamiltonian near the 
m:m + 1 resonance takes the following form: 


H = — 


GM+mi GM*to 2 GTO1TO2 


2 CL\ 


2(12 


02 


/m+l,27(ctres)^l COS \6 7X7 1 ] -T [/m+1,31 (a r es) ^m,l2ct re s]c2 COS [6 7772] 


(A.l) 


where 6 = (to + 1 )A 2 — toAi and a,;, ej, A^ and voi are the se mimajo r axes, eccentricities, mean longitudes, and longitudes of 
periastron of the two planets. The quantities /,.27 and ,/j .3 i are (Murray & Dermott[l999 p. 539-556) 


fj, 27 (a) 
/i,3l(aO 


Aj(a) 


1 . „ d \ , , , 
2 ( - 21 - a d^ )A ‘' a) 

2 ( — 1 + 2 / + 

1 r 2n cos (/</>) 
n Jo \/l — 2 a cos 4> + a 2 


(A.2) 


After expansion about the resonance center, these coefficients are evaluated at a = 01/02 = a res , where a res = [m/(m+ l)] 2 / 3 . 
Treated as functions of to, these coefficients are well approximated as 


fm+ 1,27 ~ —0.8TO 

/m+1,31 ~ 0.8to (A.3) 

as shown by |Quillen] ( |2011 ] i (and easily confirmed numerically). 

The osculating elements are not canonical. Instead we use the following as our variables: 


A i = mi\J Gm*di 
A i = Qi + tut + Mi 

Pi = mi\J Gm.*Oi(l - - ef) « A.jy 

Pi = -Wi = -Hj - uii (A.4) 


At first order in the inclinations the Hamiltonian is independent of O, the longitude of ascending node, and I. This implies 
that the inclinations need to be large enough for I 2 terms to be important for non-coplanarity to have an effect on the orbits. 
Therefore, the nearly coplanar regime is also well described by our formulation even though inclinations do not appear in the 
Hamiltonian. 

Instead of the polar set (P,p), we use the cartesian one defined as 


Xi = \j2PiCOSPi 
yi = y/2Pi sin pi 


which is more appropriate in the limit of low eccentricities. Then the perturbation Hamiltonian Hi can be written as 


H,oc4 x 


/m+l,27(a r es) , /m+1,31 ( a res) 

Xl - - -X' 2 COStf- 


v/Ai 


v/A/ 


( fm+ l,27(ares) fm+ 1 ,31 (ares) 


V v/Ai 


-y 1 




y 2 sin 6 


(A.5) 


(A.6) 


where for simplicity f m+ i 31 (a res ) = [/m+ 1,31 (a res ) ~ ^m,i2a res ]. Moreover, after expanding about the resonance center, all 
A i are held fixed at their resonant values A i res . As messy as this looks, then, the coefficients of the cos 0 and sin 6 terms are 
simply linear combinations of 2 + and y,. With an appropriate normalization, these combinations are a canonical rotation of the 
original cartesian eccentricity variables. In fact, the coefficient of the cos 6 term is a canonical momentum to a coordinate equal 
to the coefficient of the sin 6 term. This is the key - this rotation takes the 2 degrees of freedom associated with the eccentricity 
of each planets and reduces it to one degree of freedom (the conserved quantity generated here is T 2 and defined below). 
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Furthermore, since each A* only appears as the combination 9, there is a further reduction in the number of degrees of freedom 
and an associated conserved quantity which we refer to as Oi. The action associated with 6 is denoted as 0, and, 

m 

0i =-—A 2 + Ai 

m + 1 

0 = A 2 /(m + 1) (A.7) 

(see e.g. |Deck et al.|20l3l for details of these canonical transformations). We then choose units, such that actions are measured in 
terms of 0i, the Hamiltonian in terms of 0i/n 01 , where 


n 01 = 



(A.8) 


and such that time is measured in units of l/n 01 . That implies that derivatives of the Hamiltonian, which result in equations of 
motion for the variables, are time derivatives with respect to the unitless time t = in 01 . Hats will denote this set of variables. 

We define 


Xi = Xi/ \/0i 

Vi = Vi/V ®1 

§ — 1 fm+ 1,27 

A 2 ,res \/A ] .res 
1 = _ 1 /rra+1,31 

A' 2 .res V/A^^.res 

6 = J% + 5* 

5ixi + S 2 x 2 
n ~~ ~ ^ 

SiVi + hyi 
51 =-1- 

$= \{rl + s\) 

tan ib = — 
r 1 


(A.9) 


where coefficients of the rotation alluded to above can be read off the Hamiltonian and are given by 5\ and 62 - The rotated 
combinations are ri (“momentum”) and ,s-| (“position”). The momentum $ and angle if) reflect a polar canonical transformation 
of the cartesian variables r\ and Si. ( I> is related to the eccentricity as 


$ 

7 2 

C 

R 


m + 1 


a r 


m a res + C 2 (^ 2 

1 Jl 


/C^res ) 

{R 2 ei + — 2Reie 2 cos Aw) 

mi 


m 2 

| /m+1,27 (ctres) | 
f'm+ l,3l( a res) 


(A. 10) 


After the rotation and change to the polar variables $ and ni has been performed, only a single combination of the remaining 
angles 'ip and 9 appears in the Hamiltonian. A final canonical transformation to this angle and a its canonical action (equal to $) 
is performed. The last canonical transformation yields a single degree of freedom Hamiltoniarj^jof the form 


H = — r) 2 — ei(5v/2$cos^, 


(A. 11) 


7 In the |Deck et al.||2013| paper, the following Hamiltonian is defined as 


K and given in Equation 26. 
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where 


o + 4> 

3 m 2 {a res + C ) 5 

C Otres 

mi 

M~ 

$ - 50 

50 = 0 — 0 

0 _ Q-res 

m(a res + C) 



(A. 12) 


The quantity T is conserved (it is the conserved quantity found after realizing that only the combination 0 + Vj appears, and not 


each individually; refer to e.g. Deck et al. 2013). When 50 = 0, or 0 = 0, the system is at the exact period commensurability 


and a = a res . Finally, we turn now to the conserved quantity denoted 2 , found after performing the rotation. It can be written 
in terms of the Poincare variables as 


^2 


ri 


s 2 


tan 02 



S 2 x 1 - 5ix 2 


6 

&2V\ - Alj/2 

5 

£2 

T2 


(A. 13) 


T'o is conserved because the angle ip 2 does not appear in the original Hamiltonian. Note that if we wish to extend this analysis 
to second order in the eccentricities of the planets, the angle ip 2 appears in the Hamiltonian and 'l >2 is no longer conserved. These 
new “eccentricity” vector components (r,. s, ) are obtained from a linear transformation of the original (x t . yf), and <I> + T 2 = 
(Pi + P 2 )/0i since the trans forma tion is a rotation that preserves length. 

At this stage (Hamiltonian ( |A. 11 1 >), our variables are the momentum (action) $ and the conjugate angle 0. One final rescaling 
is performed to put the Hamiltonian into the form given in Equation |2]). Primes will refer to scaled quantities. We divide actions 
by the unitless quantity Q. We rescale the Hamiltonian and time t using a second parameter a, as H' = H/a and t' = ta/Q. 
Choosing a = Q 2 \/3\, and Q = (ei5/|/3|) 2 / 3 leaves us with a single free parameter (T')- The parameter Q is 


q = 4' 3 c( 


( /m+ 1,31 ^ 

1/3 5/6 

1 '-*res j 

/ + (, \JOi res \ 

\9(m+ 1)(1 + C) 2 / 

1 m 1 

\ ( Qtres C ) 5 / 


1/3 


where we have changed from to the total mass of the planets e p - e-\ + e 2 . 
The final Hamiltonian is 


(A. 14) 


h’ = - r ') 2 


V 2 $' cos 0. 


where the “time” derivatives of the variables and 0 are with respect to t' = triQ iQ|/3|. 
The conservative system is governed by the equations 


d-T 

dt’ 

dcf> 

dt’ 

dT 

dt' 


dH’ 
90 
dH’ _ 
1 )¥ ~ 


— V 2 &' sin 0 


($' _ r) 


1 


COS0 


C 


= 0 


where 


indicates the conservative evolution. 


C 


(A. 15) 


(A. 16) 


A-2. Non-conservative forces 



















16 


Deck et al. 


We now add to the conservative system the effects of eccentricity damping and migration. We parametrize these effects as 


1 dei _ 1 1 

e-% dt' Q\f3\ne 1 r ei 
1 da,i _ 1 / 2pef 1 \ 

at dt' Q\/3\n &1 V T e . r a .) 


(A. 17) 


We now convert these time variations in semimajor axis and eccentricity into time variations of <f>' and I '. First, 

1 dAi 1 ddi 1 / pef 1 \ 

A i dt' 2 cii dt' Q\/3\ne 1 \ T ei + 2 r ai ) 

1 dPi 1 dAi ^ 1 dei 1 / pef + 2 1 \ 

Pi dt' Aj dt' e, dt' Q\/3\nQ 1 \ T ei 2r ai J 


At this point, we make the assumption that \D 2 = 0 and remains zero during the evolution of the system around and in resonance. 
This greatly simplifies the equations for <f>' and (j> as we will show immediately below. That \F 2 is small when the resonance is 
encountered is consistent with the fact that the orbits are undergoing eccentricity damping and therefore have nearly circular 
orbits prior to resonance encounter. As long as the eccentricities remain low, the motion is well described by the Hamiltonian at 
O(e) used here, which is independent of A‘ 2 , and which therefore conserves \D 2 = 0. However, if the eccentricities grow too large 
as the system evolves in the resonance (since capture into resonance excites eccentricities), the 0(c 2 ) terms in the Hamiltonian 
will become im portan t, and '17 will no longer be conserved at zero. We will discuss when our assumption of \& 2 = 0 V t breaks 
down in Section [ATil 

Now, >1) + \P 2 = (Pi + P 2 )/©i, and 1) T 2 ~ 0, and 2) it remains so (i.e. the time derivative is small). Therefore, we can write 
$ = (Pi + P 2 )/0i such that 


d& 


. 1 

C]J\ +dP 2 

Z dQl ) 

dt’ 

dis 

’ 0iQ^ 

v dt' 

* dt' J 


1 

\ Pl ( 

'-pel- 2 

1 > 

+ P2 [ 

'-pe |-2 

1 ^ 

Q 2 \P\n<e> 1 

[©I \ 

Te, 

2r al J 

+ 01 l 

Te 2 

2 T a 2 ) 


$ f rn A 2 / ~2 pe\ _1_\ + A\ ( ~ 2 P e l _(*_\\ 

©1 + 1 2 V T e 2 T a2 ) 2\ T ei T ai ) ) 


(A. 19) 


Moreover, \& 2 = 0 implies that <5|Pi = A(P 2 , and therefore from $ = (Pi + P 2 )/0i and P, = A,e 2 /2 we can write 


0i 

P* 

0i 


<1> 


!+(fc 


= 7i<5$' 


$ 


1+ T 1 


<?2 


= 72 Q& 


e\ = 2?7iQ7i$' 
e 2 = 2 Q? 7 2 7 2 T>' 


0 i 

Vl ~~Ai = 

0i 

^A-r 

so the equation for the slow evolution of T' reduces to: 


m 


m + 1 (y/ctres 
m 


m + 1 


■ (j \J a r 


dt' 


dis Q\/3\n Sl 
1 


(Co + C^') 


C 0 = -7i 


2 T n 


A, + 


2rj\Ta 


■ 72 


2 T„ 


c 1 = 2 pQ ( 71 ^ ~ + 72 (to/(to + 1) ~72t? 2 ) \ 

1 ^ V T ei T e2 ) 


m/ (to + 1) 
2r]2T a2 


(A.20) 


(A.21) 


Provided \l/ 2 = 0, the angle 0 = (ra + 1)A2 -A] + A has no time evolution due to non-conservative forces. The mean 
longitudes are independent variables that are unaffected by the slow evolution of a and e. However, the angle ip docs depend 
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on di and e* explicitly, through its dependence on x t and yi. First, = 0 implies that pi — P 2 = 7r, tan^ = tanpi and 
tani/j = tan (p 2 + 7 r) = tanp 2 (recall <5i is negative while 62 is positive). From the definition ip = arctansi/ 7 ' 1 , we can derive 


d( tan ip) 


d(s 1 ) d(n) 

-tan^- 

r 1 r 1 

Sid(yi) + S 2 d{y 2 ) 


Si 

n 


r 1 


— tan 7 /) 


JicZ(ari) + S 2 d(x2) 


r 1 


2-1 


d(yi) — tan , 4>d(xi) 


(A.22) 


but d{yi)/yi = d{xi)/xi, so the quantity in parentheses is zero. Therefore the derivative of ip is zero as well: 


^d(2/i) - tanV’d(xi)^ = d{y{) ^1 - ^anp ! ) 


%i) 



= 0 


Lastly, we turn to the evolution of the single parameter T' appearing in the conservative Hamiltonian. 


dT' 

dt' 


dis 


Ao 

^t 


J-($' - se') = — - - — ^— 

dt' K dt' Qdt'\m +1 Oi 

d$' 1 d /A 2 \ 

dt' Q(m + 1) dt' \@i/ 

wFT ( a , + (A + a)*' + <?,*“) 
_I_ (± 

2 Q{m + 1)772 \r a2 


772(tn+ 1) 


- 1 


+ 




2 p f f/ 2 _ f 771 

TO + 1 \T e2 [?? 2 (to + 1) 


7i 

772 T ei 


(A.23) 


(A.24) 


A-3. Fixed point of the non-conservative system 

In total, the full equations of motion are the sum of the conservative and non-conservative pieces: 

r/cT> / ,_ (f/ 

( Ei )_ = _^ sta ^__ (Co + Cl *-) 

(E 2 )+=-(*'-n-+ r ^ 

(E3) w = ( A ° +{Al+ Co) *' +c ‘*' 2 ) <A - 25) 


We note that Q ex cj 3 and therefore Ci<J>' 2 oc e 4 /(e 2 / 3 r e ), (Co + Ai)& oc e 2 /(e 2 / 3 r e ) and A 0 oc l/(A'e 2 / 3 r e ), where 
K = T a IT e . Therefore, we can neglect the C\ term in Equation (El) above provided that e 2 <C 1, and as long as 1/AT A- e 4 , the 
Ci term in Equation (E3) is also much smaller than the Ao term. From here on, then, we ignore the Ci contribution. 

Setting all of these equations equal to zero allows us to solve for the equilibrium of the system. Equation (E3) implies the fixed 
point is given by 




eq 


Ap 

Ai + Co 


(A.26) 


Since <!>' = i’/Q, & eq = —QAq/{Cq + Ai) which is independent of Q since Ao oc 1 /Q. This implies that the equilibrium 
eccentricities are independent of e p , the total mass of the planets relative to the star. Very roughly, Co + Ai oc +, since r e <C r a , 
while QAo oc +. Therefore, $ oc tf :q oc T e /r a = 1/K. This is consistent with dropping the Ci term in Equations (El) and 
(E3), which required that e 2 <C 1 and 1/K <C e 4 . A faster rate of eccentricity damping, relative to the migration, results in a 
larger K, and a smaller equilibrium eccentricity. 

From Equation (El), we can solve for sin cp eq : 


sin <p eq 


V*Kq ^ £eq 1 

—^-Co oc-— oc- 

2Q\d\ne 1 e p n &1 T e e p n &1 T a 


(A.27) 


(note Co < 0). Lastly, we need the equilibrium value V eQ . It is a good approximation that cos (p eq = — 1 because the deviation of 
<p eq from 7 r is small for slow dissipation (n© , r„ 1 in Equation (| A.27 [> ), and the difference between cos cf> eq and -1 is quadratic 
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in this deviation. In that case. Equation (E2) of Equations (|A.25|> yields 


r 


/ 

eq 


$1 


eq 


1 

V^q 


(A.28) 


In the conservative problem, cf> eq = n. Equation (E2) of Equations ( |A.25[ i does not contain dissipative terms, and therefore 
the relationship between T' eq and 4>' eq is approximately unchanged by the dissipation since cos <p eq is approximately unchanged 
as well. That implies that the fixed point of the dissipative problem (</> eg , T' eq ) lies very close to the fixed point of the 
conservative problem with ((f) eq = 7t, & eq ) with a proximity parameter of l’ 7 = r' This is the fixed point labeled X\ in Figure 

m 

Lastly, so as not to introduce extra parameters in the main text, we wrote the dissipative terms appearing in the equations for 
<!>' and T' as 


d-f-' 

dt’ 

dF 

dt' 


= C, V 


dis 


ao 

T a 


Co pa i 

T~e At,e 




(A.29) 


dis ' a 

where the parameter ai here has nothing to do with the semimajor axes. By comparison with Equations ( | A.25 [ > one finds that 

Cq T e 


Co = 

a o = 

a\ = 


Q\/3\n ei 

A 0 T a 

Q\(d\ne 1 

AiT ate /p 

Q|/3|n 01 


(A.30) 


where we had anticipated that the C\ term is negligible. 


A-4. Linear stability analysis of the fixed point 
The matrix whose eigenvalues we wish to determine is 



-d(<^)/dV 


r)/dfi 

oJB- 
& o 

_i 

M = 

d($)/d*' 

»(| 


d(^)/dT' 


[ d^/dV 

9 (f t 

t)M 

d(%r)/d rj 

evaluated at the equilibrium. We write this symbolically as 






kA 

B O' 



M = 

C 

nD 1 




nE 

O 

O 



where k denotes a term proportional to 1 / (n@, r„) coming from the non-conservative addition to the equations of motion and 

= + = 2Q|/3|n 01 

B = - y^cos fi eq « ^2$^ + 0{k 2 ) 

C = - 1+ (2i^ Cm<P ‘- *~( 1+ (2*6^) + ° { * 2) 

A 1 . , Co 

= Ql/SIne, {Al + Co) <A ' 31> 

We will solve for the eigenvalues perturbatively in n. At this stage, we have taken advantage of the fact that the correction to 
cos (f> eq = — 1 is second order in n. 

The equation for the eigenvalues is 

0 ={kA - A ){kD - A)(—A) + BnE + A BC 


(A.32) 
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We seek solutions of the form A = (nao, nai ± i(a 2 + ko^)). Plugging in the first root yields an equation for a 0 . Plugging 
in the second and third roots yield 4 equations, for the real part and imaginary part of each root. Only two of these are distinct. 
Truncating at zeroth order in k yields an equation for ct 2 , while at first order we find two equations for a\ and 0 - 3 . We find that 


ot 0 


Ctl 


«3 


E 

A + D ( E 
2 + 2 C 

-BC > 0 = w 2 

0 


(A.33) 


where u = \J\BC\ is the unperturbed frequency (of the conservative system) evaluated at the modified fixed point ($' = 

Kq Aeq ~ 7t)- 

We assume for simplicity that m■ ~ (to + 1), R ~ 1 and a res ~ 1 (the compact approximation mentioned in the main text). 
Then: 


71 = 

Vi = 

72 = 

V 2 = 1 + C 


1 


C + i 



C 


C + i 


Q = 
$ = 


r 2/3 C 

p (i + C ) 2 
C 2 
2(1 + 0 



1/3 


a 2 = el + e\ — 2eie 2 cos Aru 


(where we have used Equation (|A. 10|>), so that 
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^0 

+1 


2 1 1 — C 1 

C + 1 T e 2(1 + 0 T a 



2e 2/3 T a 

2 p 1 


TO(1 + C) 2 T a ,e 


where 


1 

1 

1 

— 

= - 
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T a 

T a , 2 

T a , 1 

1 

1 

c 

— 

= - 

+ — 

T e 

Te, 1 

T e , 2 

1 

1 

( 2 a 

7~a,e 

T e ,l 

R 2 T e , 2 


(A.34) 


(A.35) 


(A.36) 


Anticipating the discussion regarding r a e , we do not make the compact approximation only for this parameter. 

The equilibrium value of & eq can then be evaluated and used to determine the equilibrium value of a eq by referring to Equation 
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( |A.34[ ), as 


& = 
eq 

0 

\ 0 

Ai + Co 

4m T a 
C+l r e 

_2 

a eq ~ 
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2m r a | 2 p 

T a 


9m" 

O.S 2 e 2 


1/3 


4 p 


= <y\ 


c 


eg 


(1+C) 2 


Q 2(1 + C) 2 


(A.37) 


where we have assumed that r a /§> r e and hence ignored the second term in Co- To obtain the equilibrium eccentricity without 
coupling between the semimajor axis and eccentricity evolution, set p = 0. Also notice that if r e/i —>• oo (the limit of no 
eccentricity damping), both r a e and r e diverge. In this case there is no equilibrium eccentricity. 

The stability of the fixed point is determined by the sign of a t and the sign of tt (l . First, the sign of a ( , s[cti], is given by: 


a[ai] = s 
0 >«- 
Ai + Cq ~ — 


C' 0 (1 + (2$')- 3 / 2 )-(A 1 + C7o) 


C 0 ( 2 $' )" 3 / 2 - A x 


2 1 
C + 1 Te 

2 p 1 2 1 

m(l+Q 2 T a ,e C+lfe 

where we have assumed r e T a (K /g> 1). Then 


a[ai] = s 


2 1 / 2 t„ 2 p 


C+lr e VC + l T e rn{ 1 + C) 2 


3 ^ 2 0.8 y/me p 2 p 1 
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Rearranging, a\ > 0 when: 


3 pm r e 


£p < 623/2 Ta e 


(1 + C) 2 


i(C +1) +p^, 


I - 

3 / 2 V 7V, 


3/2 


(A.38) 


(A.39) 


(A.40) 


Here 6 = 0.8 m. We remind the reader that the criterion given here has employed the compact approximation. If more accuracy 
is needed (likely only for the 2:1 resonance, refer to f igure hj, return to the full expressions for A, D, C and E, employing the 
full expressions for 71 , 72 , Vi , and Q, to e valuate o i (or whatever else) directly. 

Importantly, the criterion Equation ( |A.40| ) can never be satisfied if T a _ f . < 0, or if 

Te,2 < C 2 aT e ,i/f? 2 

r e,2 < C 2 Te,l- (A.41) 

(Again, the factors of R and a are only important for the 2:1 commensurability). 

In the case where p = 0, 


s[ai] = s 


C + l r e\C+l T ! 


3/2 


0.8 y/me 


= S 



. r e \T e ) _ 


(A.42) 


which is always negative for convergent migration T a > 0 and eccentricity damping T e > 0. This shows the stability of the fixed 
point when the eccentricity-semimajor axis coupling term is not present. 

Next, we turn to the sign of ao- Since C < 0, s[ao] is determined by s[-E]: 


s[a 0 ] — s[.E] — s[A_! + Cq] ~ s[— 


P 


w(l + C) T a>( 


- 1 


(A.43) 


Note, however, that for & ea to be a positive quantity (for the fixed point to exist), its denominator must be positive. From the 
definition given in Equation ( |A.37| ), we see that this means that (canceling common positive factors between the two terms in the 
denominator) 1 + [p/m){r e /T a ^ e )( 1 4 - £) -1 > 0- In that case, s[«o] < 0 always. 
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A-5. Limiting case of a test particle 
A-5.1. Inner test particle 

We consider the case where the inner planet is a test particle, moving outwards, towards a massive planet on a fixed circular 
orbit. Then ^ —>• 0, r a> 1 —> — r a ,i, t 0i2 —» oo, r e , 2 —► 00 and 


The criterion is 


1 _ 1 

Aj. t a l 

1 _ 1 _ 1 

7"a,e As Ts,! 


£2 < 


3 pm 1 

£>2 3 / 2 (to + p) 3 / 2 



(A.44) 


(A.45) 


Plugging in p 


1 and converting from r 0( i to r n i yields: 


£2 < 


BV3(1 + m) 3 / 2 ) 



(A.46) 


The equilibrium value for eccentricity is given by: 

2 1 Te,l 1 Te,l 

1 2(p + m) T ajl 3(p + to) t„ ! 

The equilibrium values of sin (f> is given as: 


1 T e , 1 

J= i 3(1 + to) l 


(A.47) 


sin = 




Cn 


69 2Q|/3|n 01 

Q|/3]n 01 = n 2 £p /3 (0.8 2 m 4 3) 1/3 


sin (f> eq = - 


ei 


T ei iniepS 


(A.48) 


which was found by Goldreich & Schlichting (2014|>. 


In this case, r Q l —► 00 , r 6j 1 —> 00 , and C 


A-5.2. Outer test particle 
00 . The equilibrium eccentricity is given by 


Te, 2 


1 


2(?n — p) 2K (to — p) 


(A.49) 


We emphasize that this estimate of the equilibrium eccentricity is an approximation. It does appear that for some parameter 
choice s (to > p), the equilibrium eccentricity does not exist. This means that there is no fixed point in resonance. According 
to |Tanaka & Ward| ( 2004) l, p < 1, so since to > 1 this issue may not arise in practice. However, in the case of p = 1, we have 
confirmed that the full expression for the equilibrium eccentricity does not yield a diverging equilibrium for to = 1 in the case of 
the outer test particle (see Figure|7]i. 

The criterion for over stability is 


£i < — 


3 pm 


_ 1 ( T e,2 

B2 3 / 2 (to — p) 3 / 2 \ r a> 2 


3/2 


(A.50) 


This implies that a \ is never positive. Note that we assume to > p to get the equilibrium, which also implies that o: 0 < 0, and 
that the criterion for over stability ei cr u is a real number. 


A-6. The assumption that 4r 2 = 0 

Since we have prescribed eccentricity damping, the assumption that T 2 is zero when the resonance is reached is a good one. 
However, if the system is captured into resonance, the eccentricities grow as 4> grows, and the condition that T 2 = 0 requires that 
this growth preserves the orientation of the orbits and the ratio of the eccentricities. This is true because when 4r 2 = 0, r 2 = 0 
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Fig. A. 1.— Comparison of the predicted equilibrium eccentricities in the case of ^2 = 0 (solid) with observed equilibrium eccentricities (dashed) for the 2:1 
mean motion resonance for different values of the planetary mass ratio £ (as labeled). The magnitude of the eccentricity damping determines where the system 
halts along the dashed curve. 


and s 2 = 0, which in turn imply that 


ZU\ — ZU 2 

Pi 

or 

e 2 = ei C,y/a/R 

tan-)/’ = tanpi (A.51) 



with R = |/ 27 /(/ 3 i - 2a<S m ,i)|. 

It cannot be true that \I/ 2 = 0 genetically for large eccentricities. For example, it has been demonstrated numerically and 
analytically that for the 2:1 resonance the center of resonance changes as the eccentricities grow, from anti-alignment of the 
orbits (v o\ =■ m 2 + 7 t) to exact alignment = n7 2 ) or asymmetric libration, depending on the mass ratio of the planets to each 
other ( Lee & Peale|2002] [Beauge et al. 2003[ >. This change implies \k 2 deviates from zero. 

To test how large the eccentricities can be before the approximation breaks down, we performed numerical tests, driving 
systems with T' 2 (f = 0) = 0 into resonance without eccentricity damping. This leads to an unchecked increase in <!>' and in 
<7 k, 'J+ e\ — 2eie 2 cos Aw. If this growth preserves \J/ 2 = 0, the eccentricities will grow along a track of constant e\ /e 2 . 
Eccentricity dam ping, if included, would halt the system at a particular value of a along this track. 

In FigurelA.il we show, as a function of £, the comparison between the analytic prediction with = 0 and the results of 
the numerical integrations for e\ and e 2 as the system is driven into the 2:1 mean motion resonance. We show the same for the 


3:2 resonance in Figure [A2| 
assumption that T' 2 k, 0 to 1 


How low eccentricities need to be (or equivalently, how effective the damping must be) for the 
Be good is also a function of the resonance integer to. Especially for the 2:1 resonance, it is clear that 
the true evolution drives \k 2 away from zero as eccentricities grow, but for the 3:2 resonance the prediction based on T 2 = 0 is 
very good across the entire range studied (ei, e 2 ) £ [0, 0.1]. For the 2:1 resonance in particular, depending on the equilibrium 
eccentricities, the predictions of the theory may be only roughly correct, and they may be wrong entirely in the case where the 
resonance center corresponds to alignment or asymmetric alignment of the orbits. 
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Fig. A. 2.— Comparison of the predicted equilibrium eccentricities in the case of '$'2 = 0 (solid) with observed equilibrium eccentricities (dashed) for the 3:2 
mean motion resonance for different values of the planetary mass ratio £ (as labeled). The magnitude of the eccentricity damping determines where the system 
halts along the dashed curve. 



